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ABSTRACT 


Transient  analysis  deals  with  the  problem  of  numerically 
determining  the  poles  si»s2»***  associated  with  the  Laplace  trans¬ 
form  of  a  given  real  transient  y(t)  «  £  av  exp(svt)  (where 
Re  sv  <  0  for  each  v) ,  from  a  knowledge  of  its  samples 
y-j.  =  y(k  At),  k  =  0,1,...  .  In  practice,  the  samples  y0,y1,y2>... 
are  usually  contaminated  with  noise,  and  this  serves  to  limit  the 
effectiveness  of  the  computational  schemes  of  Prony,  Bellman,  Jain, 
Van  Blaricum,  etc.  which  have  been  developed  for  extracting  the 
first  few  sv's.  The  performance  of  these  algorithms  can  be  greatly 
enhanced  if  the  data  yQ,y^,y2»...  is  first  subjected  to  a  suitable 
sequence-to-sequence  transformation.  Any  such  linear  pole 
preserving  transformation  must  have  the  simple  form 

wk  =  Wk  +  Ylyk+1  +  y2yk+2  +  k  =  where 

Y(z)  -  Yq  f  Y^z  +  "i +  •••  is  analytic  and  zero  free  in  the  unit 

disc  |z|  <1.  In  most  cases  of  interest,  y(z)  may  be  chosen  so 

as  to  greatly  suppress  the  effects  of  noise. 
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PREPROCESSING  TECHNIQUES  IN  TRANSIENT  ANALYSIS 
David  V.  Kammler 

ABSTRACT 

We  consider  the  problem  of  finding  the  poles  z^,  z^,..., 

z  associated  with  the  z-  transform  of  the  sequence 
n 

y=  (y0>y1>y2* • • )  of  samples 

yk  =  Y(kT),  k  =  0,1,2,... 

of  the  transient 
n 

Y(t)  =  avexp(svt),  t  >  0 
v=l 

emitted  by  a  given  n-th  order  linear  system.  In  principle, 
exactly  the  same  poles  can  be  extracted  from  the  sequence 
u  =  T  y  when  T  is  a  sequence-to-sequence  transformation  of 
the  form  r  =  V(E)  where  E  is  the  shift  operator  and  Y(z)  is 
analytic  and  zero  free  on  the  unit  disc  |z|  <  1  .  Such  a 
preprocessing  operator  T  can  be  chosen  so  as  to  suppress 
additive  noise  or  to  selectively  enhance  one  or  more  of  the 
poles  without  annihilating  the  others.  Using  such  prepro¬ 
cessing  operators  we  obtain  a  common  conceptual  framework  for 
all  of  the  previously  used  schemes  for  transient  analysis 
(including  those  of  Prony,  Van  Blaricum  &  Mittra,  and  Jain) 
and  we  provide  a  theoretical  basis  for  several  promising  new 
algorithms. 
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1.  INTRODUCTION  TO  TRANSIENT  ANALYSIS 

Let 

(1)  Y(t)  =  £  aveSvt’  t  >  0 

v=l 

denote  an  n-th  order  approximation  to  the  transient  emitted 
by  a  physical  system  in  response  to  some  initial  excitation 
c  f ’ .  [  1  ] .  (In  the  appendix  we  present  a  mathematical  model  for 
a  damped  vibrating  string  which  serves  to  illustrate  the  kind 
of  phenomena  we  wish  to  analyze.)  Ve  consider  the  problem 
of  numerically  determining  n  and  the  complex  frequencies 
Sl,"'',Sn  ^rom  a  knowledge  of  the  sequence  y=(y^, y^ , y2 , . . . ) 
of  uniformly  spaced  samples 

(2)  yR  =  Y(kT),  k  =  0,1,2,... 

where  T  >  0  is  the  sampling  interval.  We  shall  assume  that 
Y  is  real  valued  so  that  the  av’s  and  sv's  occur  in  complex 
conjugate  pairs,  and  ve  assume  that  I  is  nondegenerate  in 
the  sense  that  the  s  '  s  are  distinct  and  the  av  's  are 
nonzero.  We  further  assume  tnat  the  sv's  all  lie  in  the 
left  half  plane  so  that  Y(t)— >  0  as  t  — >  3C  .  By  substituting 
t=kT  in  (1)  we  see  that  the  samples  have  the  representation 

(3)  yk  =  22  avzv  ’  k=0,l,2,... 

v=l 

where 

,  ,  svT 

(4)  zv  =  e  ,  v=i, .  .  .  ,  n 
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with 


(5)  lz^i  ^  1* 

Since  s^,...,sn  occur  in  complex  conjugate  pairs  the  same  is 

true  of  z.,...,z  so  that  the  coefficients  c_,c,,...,c  of 
In  O  1  n 

the  corresponding  characteristic  polynomial 

(6)  cU)  =  C()  +  +  ...  +  cnz'1,cn(2-21)...(z-Zn),cn€K,cn^0, 

are  all  real  valued.  Ve  shall  assume  that  the  sampling 
interval  T  is  sufficiently  small  so  that  z^,...*z  are 
distinct  and  so  that  the  inversion  of  (4) 

(7)  sv  =  T  1  fnzv,  v=l,...,n 

can  be  effected  unambiguously  without  undue  regard  to  aliasing. 

Ve  shall  find  it  convenient  to  introduce  the  shift 
operator  E  which  is  defined  so  that 

(8)  E(y0,y1,y2, . . . )  =  (y : , y2 > y 3 > . . . ) 

or  equivalently 

Eyk  =  yk+l’  k=0’1»2» •  •  •  » 

and  we  shall  let  I  denote  the  corresponding  identity 
operator.  Ve  observe  that 

2  2 
E( 1 , z, z  ,...)  =  z  ( 1 , z, z  >...) 
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>  i.  ■  ■ 


so  that  the  application  of  the  operator  E-zl  annihilates  the 
power  sequence  (l,z,z  ,...).  This  being  the  case,  since  (3) 
expresses  y  as  a  linear  combination  of  n  such  power  sequences 
corresponding  to  z^,...,z  we  must  have 

(9)  (E-z1)...(E-zn)y  =  O. 

By  using  (6)  we  may  rewrite  (9)  in  the  form 

c(E)y=(c()I+c1E+.  .  . +cnEn)y=cQy+c1  (Ey)  +  ..  .  +cn(Eny)=0 
or  equivalently 

(10)  iJ  c  =  O 

where  the  data  matrix 
y0  yl  yn 

(11)  y  =  yl  y2  ‘ • •  yn+l 

•  •  • 

•  •  • 

•  •  • 

J  n  X 

has  the  columns  y,  Ey,  . . . ,  E  y  and  where  c  =  ( cq> c ^ > • • • > cn )  . 

Since  we  have  assumed  that  t?  '  av ' s  are  nonzero  and  the  zv's 
are  distinct,  the  first  n  columns  of  are  linearly  indepen¬ 
dent,  d  has  rank  n,  and  (apart  from  a  scale  factor)  the  null 
vector  c  is  uniquely  determined., 

These  observations  suggest  the  following  approach  to  the 
transient  analysis  problem.  Given  the  samples 
y  =  (yQ,y1,y2f . . . )  we  attempt  to  make  some  slight  over¬ 
estimate  of  the  system  order  n  and  form  the  data  matrix  (ll). 
We  numerically  investigate  the  null  space  of  3  and,  if 


necessary,  reduce  n  so  as  to  obtain  an  essentially  unique 
null  vector  c.  After  effecting  the  factorization  (6)  and 
using  (7)  we  then  obtain  the  desired  s^'s.  Almost  all  of 
the  presently  used  noniterative  schemes  for  computing  the 
sv's  from  the  y^ 1 s  fall  within  this  conceptual  framework. 
For  example,  when  we  use  the  well  known  method  of  Prony, 
we  ignore  all  but  the  first  n  rows  of  ij  and  after  arbi¬ 
trarily  setting  cfi=l  we  solve  the  resulting  system  of 
linear  equations 


y0  yl 

*  *  *  yn 

co 

0  ' 

yl  y2 

• 

*  *  *  yn+l 

• 

• 

0 

• 

• 

• 

c  i 

n-l 

• 

• 

yn-l  yn 

y2n-l 

1 

0 

to  obtain  the 
vector  c. 

Perhaps 
vector  for 
noise)  is  to 
length  of  yc 
length  of  c, 
quotient 

T  T 

(13)  c  aac 
T 

c  c 


remaining  components  of  the  approximate  null 

the  most  natural  approach  to  finding  a  null 
(especially  when  ij  has  been  contaminated  with 
determine  c  so  as  to  minimize  the  Euclidean 
subject  to  some  normalization  of  the  Euclidean 
or  equivalently,  to  minimize  the  Rayleigh 


(In  so  doing  we  avoid  changing  an  intrinsically  homogeneous 
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A 


problem  into  an  inhomogeneous  one  as  in  (12)  by  the 
imposition  of  a  constraint  cn=l  on  one  of  the  components  of 
c. )  Any  minimizing  o  for  ( 1 3 )  must  be  an  eigenvector 
belonging  to  the  smallest  eigenvalue  of  the  symmetric 
nonnegative  semidefinite  (n+l)  x  (n+l)  matrix  ^  y , 
cf . [8 j p. 266  ]  .  In  the  ideal  noise  free  setting  where  the 
system  order  n  is  known,  both  H  and  y  will  have  rank  n, 

the  minimum  eigenvalue  of  if  y  will  be  zero,  and  the  corres¬ 
ponding  eigenvector  c  (unique  to  within  a  scale  factor)  will 
be  the  desired  null  vector  of  3  .  In  practice,  we  do  not 
know  the  system  order  and  we  can  only  approximate  the  elements 


(yTs>io 


-  £ 

k=0 


=  0,1,. 


due  to  noise  contamination  of  the  data  and  to  the  use  of 

finite  precision  arithmetic  on  our  computer  so  our  computed 
T 

y  y  will  not  have  either  a  zero  eigenvalue  or  a  null  vector. 
Nevertheless,  by  using  well  known  numerical  methods  [7]  we  can 
effect  the  eigenvalue  decomposition 


(15)  y  y  = 


T  ,  T 
ro  +  Xl',lTl 


.  .  +  X  Y  Y 
n  n  n 


where  Xq  >  >  .  .  .  >  X  >  0  are  the  eigenvalues  and 

V0,Vl,'*",Vn  a  corresPoriding  orthonormal  set  of  eigenvectors 
T 

for  b  S  •  If  we  purposely  use  a  value  of  n  which  is  a  bit 
too  large,  we  can  analyze  the  distribution  of  the  X^ 1 s  and 
thereby  ascertain  the  correct  system  order,  e.g.,  when  n  is 
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correctly  chosen  we  might  expect  to  have  ^  >  >  X  and 
Xn  w  0.  Tne  corresponding  is  then  a  good  choice  for  e. 

This  approach  has  been  successfully  developed  by  VanBlaricum 
and  his  coworkers,  [ 1 , 9 ] - 

As  an  alternative  to  the  eigenvalue  decomposition  (15), 
we  can  use  a  singular  value  decomposition  of  $  to  obtain 
the  c  which  minimizes  (13).  The  data  matrix  i)  can  be 
represented  in  the  form 

U6>  if  =  OoV'c?  +  °iu1tiT  +  •••  +  <Wn 

where  Oq  >  cr^  >  .  .  .  >  on  >  0  are  the  singular  values  of  , 

where  v„,  ,  ...»  ▼  are  orthonormal  vectors  from  R.n+^,  and 

(J  1  n 

where  u^,  un  are  orthonormal  sequences.  (Indeed, 

Vq>  .  .  .  ,  vn  are  again  an  orthonormal  system  of  eigenvectors 

of  as  in  (15)  with  k=0,l,...,  n  and 

for  each  k  for  which  0.)  Tne  and  Vj^'s 

can  be  computed  directly  from  the  elements  of  the  data 
matrix  without  first  f orming  y ^y,  cf.[2,3].  The  distri¬ 
bution  of  the  a ^ ’ s  can  be  used  in  the  same  way  as  that  of 
the  X^ ' s  to  help  determine  the  system  order,  and  the  right 
singular  vector  serves  as  an  approximate  null  vector  c 
f or  d  . 

In  principle,  the  singular  value  decomposition  (16) 
should  result  in  a  somewhat  better  conditioned  estimate  of 
the  approximate  null  vector  c  than  the  eigenvalue  decomposi¬ 
tion  (15).  Assuming  a  .  >  a  ,  it  can  be  shown  that  when  M 

6  n-1  n  ^ 
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is  replaced  by  y  +  £  in  (16)  the  resulting  perturbation  Ac 
which  results  in  the  normalized  approximate  null  vector 


c  =  t  satisfies  the  bound 
n 

(17)  11  h  <  °_o  .  ||£||,  4  0(11  £11,) 


T  T 

whereas  when  y  y  is  replaced  by.  y  y  +  'bJ  in  (15)  we  have 
the  corresponding  bound 


(18) 


II 

if 


aell2 
•  I!  2 


<  X0  .  IIU  II 2  +  o(HU|l2), 


with  both  bounds  being  sharp.  The  condition  number 

1  2 
(19)  h  =  0  -  a0 

L  X  ,  -  X  2  2 


associated  with  the  null  vector  computation  based  on  (15)  is 
usually  many  orders  of  magnitude  larger  than  the  corresponding 
condition  number 


(20)  ks  =  CT0 

a_  ,  -  c 
n-l  n 

associated  with  the  null  vector  computation  based  on  (16). 

On  the  other  hand,  when  using  the  eigenvalue  decomposi tion 
we  must  store  only  one  copy  of  (a  suitable  truncation  of) 
the  data  sequence  y  =  (yq> Yj » Y2 > • • • )  to  use  in  computing  the 
elements  (14)  of  j  ,  and  this  symmetric  (n+l)x(n+l) 

matrix  is  then  used  as  the  input  to  a  routine  which  performs 
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the  eigenvalue  computations.  In  contrast,  we  must  assemble 
(which  requires  n+1  times  as  much  storage  as  y)  to  use 
as  the  input  for  a  singular  value  decomposition  analysis. 
When  the  system  order  n  is  large  and  when  many  rows  of 
are  known  and  available  for  use,  this  additional  storage 
requirement  for  the  singular  value  decomposition  may  very 
well  prohibit  its  use. 

The  above  formulation  of  the  transient  analysis  problem 
has  been  strongly  influenced  by  a  recent  paper  of  Henderson 
[4]  and  by  Volume  I  of  the  technical  report  [l]  of  Auton 
and  Van  Blaricum.  An  unusually  complete  annotated  biblio¬ 
graphy  of  related  papers  and  technical  reports  is  given  in 
Volume  III  of  [ 1 ] . 


i 
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2.  PREPROCESSING  WITH  THE  TAIL  SUM  OPERATOR 


Suppose  that  we  are  given  a  transient  sequence 
y=  (y0»y1>y2»- • • )  having  components  (3)  with  ai>---»an 
being  nonzero  and  with  being  distinct  points 

within  the  unit  circle,  and  that  we  wish  to  determine  the 
zv's  numerically.  Although  it  is  possible  to  process  tne 
yk '  s'  directly  as  described  in  the  previous  section,  it  is 
often  advantageous  to  use  the  y^ ' s  to  generate  an  auxiliary 
sequence  a  =  ( u^, u^ , U2 . . . )  from  which  we  subsequently 
extract  the  zv 1 s  .  This  procedure  is  known  as  preprocessing. 
Before  giving  a  more  precise  definition,  we  shall  consider 
a  specific  example. 

Ve  define  the  tail  sum  operator  S  to  be  the  sequence- 
to-sequence  mapping  for  which 

(21)  Syk  =  J3  yk+i  *  k=0 ,1,2,...  , 

t=0 

i.e.,  to  generate  the  elements  of  n.  =  Sy  we  compute 

(22)  »k  =  yk  +  yk+1  +  yk+2  +  •••  -  *=0,1,2,... 

directly  from  the  elements  of  y.  Using  (21)  we  see  that 

Szk  =  £  zk+i  =  Zk/(  1-z) 

1=0 

so  that 

S ( 1 , z, z  ,...)  =  (l-z)  .  ( 1 , z , z  ,...)  , 
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2 

i.e.,  the  power  sequence  (l,z,z  ,...)  is  an  eigenvector  of 
S  corresponding  to  the  eigenvalue  (l-z)-*  for  every  choice 
of  |z|  <  1.  This  being  the  case,  if  y  has  the  representation 
(3)  (so  that  y  is  a  linear  combination  of  n  such  power 
sequences),  then  u=Sy  has  the  components 
n  #  , 

(23)  u^  =  ^2  av  zv  »  k=0,l,2,... 

v=l 

where 

(24)  ay  =av/(l-zv),  v=l,...,  n. 

¥e  thereby  see  that  the  sequence  u=Sy  has  exactly  the  same 

poles  z^,...,zn  as  y,  so  it  is  possible  to  extract  these 

poles  from  u  as  described  in  the  previous  section.  Using 

*  . 

(24)  we  see  that  the  ratio  av/av  is  large  when  z^  is  near  1, 

and  thus  the  pole  zv  is  more  strongly  represented  in  u  than 

in  y  when  this  is  the  case.  Moreover,  we  might  expect  the 

summation  (22)  to  suppress  some  of  the  effects  of  any  noise 

which  may  have  contaminated  the  y^'s.  This  being  the  case, 

we  might  reasonably  hope  to  extract  slightly  more  accurate 

poles  from  u  than  we  could  obtain  directly  from  y. 

Of  course,  if  a  single  application  of  S  tends  to 

suppress  the  noise  and  enhance  the  poles  which  lie  near  z=l, 

the  repeated  application  of  this  operator  might  very  well  be 

expected  to  do  an  even  better  job.  The  application  of  S  to 

2 

the  sequence  Sy  gives  the  sequence  S  y,  etc.,  with  each  of 
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the  sequences  y,  Sy,  S  y, . . .  having  exactly  the  same  poles 
zl'‘*',zn*  ^his  suggests  three  different  schemes  for 
computing  these  poles.  First  of  all,  we  might  select  some 
fixed  p=0,l,2,...  and,  following  our  earlier  approach, 
attempt  to  find  a  null  vector  o  for  the  data  matrix 


(25) 


spy0  sPyi  ...  SPyn 

sPyl  S?y2  sPyn+l 

Spy2  SPy3  ...  sPyn+2 


(which  reduces  to  (ll)  wnen  we  take  p=0  and  I.)  Upon 

effecting  the  corresponding  factorization  (6)  we  obtain  the 
zv's  .  The  second  approach  is  suggested  by  the  observation 
that  in  the  absence  of  noise  the  exact  null  vector  c  is 
orthogonal  to  the  first  row  of  the  matrix  (25)  for  each 
p=0,l,2,...  .  We  might  therefore  seek  a  null  vector  c  for 


the  matrix 


(26) 


y0  yl 

Sy0  Syi 

2  2 

s  y0  s2yi 


Sy 

q  2 

s  y 


n 


n 


and  again  obtain  the  zv ' s  from  the  resulting  factorization 
(6).  Finally,  by  repeatedly  applying  S  to  (3)  we  see  that 
n  , 

(27)  Spyk  =  ayz*  [ ( l-zy) ]p,  p=0,l,2,...,  k=0,l,2,... 

v=l 
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the 


and  thereby  observe  that  for  each  fixed  k=0,l,2,... 

sequence  yfc  ,  Sy 

vv  =  (1-.  r1,  v=l,...,n  (which  lie  in  the  half  plane 

Re  w  >  1/2  when  |zv|  <  1  for  each  v. )  This  being  the  case, 

T 

if  we  find  a  null  vector  d  =  (d^, d ^ . . . , dn )  for  the  matrix 

y  0  Sy  0  ‘  •  S  y  0 
yx  Syi  ...  Snyi 

(28)  y2  Sy2  ...  Sny2 

•  •  • 

•  •  a 

a  a  a 

and  effect  the  corresponding  factorization 

dQ  +  d^w  +  ....  +  dRwn  =  dn(w-w1)  ...  (w-vn) 

we  will  have  wv  =  (l-zv)_1  or  equivalently  zy  =  (wy-l)/wv  , 
v=l,...,n  (after  a  suitable  permutation  of  the  indices.) 

This  third  approach  is  equivalent  to  the  pencil-of-functions 
method  of  Jain  and  his  coworkers,  [5,6]. 

If  j  has  the  exact  representation  (3)  and  all  compu¬ 
tations  are  performed  without  error,  these  three  approaches 
all  yield  the  same  zv's.  In  practice,  however,  the  system 
order  may  be  infinite  (with  n  of  the  poles  being  dominant 
and  the  influence  of  the  others  being  small  but  nonzero), 
the  sampled  y^ ' s  may  be  contaminated  with  noise,  and  finite 
precision  arithmetic  is  used  to  carry  out  the  computations. 
The  accuracy  of  the  computed  zy's  thus  depends  on  which  of 


S2yk' 


has  the  associated  poles 
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/ 


the  three  matrices  (25),  (26),  (28)  we  use  and  on  whic 
method  we  use  to  determine  the  approximate  null  vector 
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3.  PREPROCESSING  OPERATORS 


Ve  shall  nov  generalize  the  results  of  the  previous 
section  so  as  to  include  preprocessing  schemes  other  than 
those  based  on  the  tail  sum  operator  S.  Let 

(29)  u  =  r  y 

be  a  linear  sequence-to-sequence  transformation.  Ve  may 
think  of  u,y  as  being  (column)  vectors  with  r  =  [  Y ^  )  being 
a  matrix  so  that 


(30)  uR  =Tyk  =  YkQy0  +  Ykly,  +  \2y2  *  ...  . 


Ve  would  like  to  impose  restrictions  on  the  y,  ' s  which  will 

K  C 

insure  tnat  V  is  defined  on  the  whole  set  of  rapidly  decaying 
transient  sequences 

(31)  J  =  {  (yQ>y1 >y2>  •  •  •  )T  :  iim  sup  I 1 1  //k  <  l  }  . 

(Vhen  y^  has  the  form  (3)  with  av^  0  and  | zv |  <  1  for  each 
V=1 , . . . , n  we  find 

lira  sup  |yR  | 1//k  =  max  {  |  Zj  | ,  .  .  .  ,  |  zn  |  }  <  1 

so  that  O'  includes  all  of  the  transient  sequences  which 
might  arise  from  any  stable  finite  order  system  we  might 
wish  to  study.)  Moreover,  to  be  useful  as  a  preprocessing 
scheme  we  must  insist  that  T  be  pole  preserving  in  the  sense 
that  zv  is  a  system  pole  associated  with  u  =  fy  if  and  only 


if  zv  is  also  a  system  pole  associated  with  y,  and  this 
serves  to  further  restrict  the  Yk^'s.  The  resulting 
sequence-to-sequence  mappings  are  cnaracterized  by  the 
following 

THEOREM  1 :  A  necessary  and  sufficient  condition  for  T 
to  be  a  linear  pole  preserving  sequenc e-to-sequenc e  mapping 
of  3"  into  CT  is  that  T  have  the  form 

(32)  r  =  Y(E)  -  Y0  +  YXE  +  Y2E2  +  ... 
i .  e.  , 

(33)  pyk  =  Y0yk  +  Y^j^  +  Y2yk+2  +  •••  »  k=0,l,2,... 
where 

(34)  Y ( z )  =  y0+Y1z  +  Y2z2  +  . . . 

is  a  zero  free  analytic  function  on  the  unit  disc  |z|  <  1  . 

Before  proving  this  result,  we  point  out  that  (32)— (33) 
imply  that  T  has  the  banded  upper  triangular  matrix 
representation 


and  since  the  radius  of  convergence,  R,  of  (34)  is  given  by 
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the  Cauchy-Hadamard  formula 


(36)  R  1=  lim  sup  |Yk|1//k  , 

the  analyticity  of  y(z)  on  |z|  <  1  is  equivalent  to  the 
requirement 

(37)  lim  sup  |Yk|1//k  <  1  . 

The  hypothesis  that  y(,z)  be  zero  free  is  not  easily  trans¬ 
lated  into  a  simple  condition  on  the  Yk ' s  . 

To  illustrate  the  theorem,  we  first  note  that  the 
tail  sum  operator  S  of  (22)  has  the  representation 

(38)  S  =  Y(E)  =  I  +  E  +  E2  +  . . . 
where 

(39)  Y(z)  =  l  +  z+z2+  ...  =  1 / ( 1  —  z ) 

is  clearly  analytic  and  zero  free  on  |z|  <  1.  Likewise, 
the  local  smoothing  scheme 

<4°)  uk  =  yk  +  yk+1  +  ...  +  yk+K.1 

which  results  from  the  operator 

(41)  SN  =  Yn(E)  =  I  +  E  +  E2+  ...  +  EN_1 
with 

(42)  Yn(z)  =  1  +  z  +  z2  =  . . .  +  zk_1  =  (l-zk)/(l-z) 
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meets  the  hypotheses  of  the  theorem.  On  the  other  hand,  the 
weighted  smoothing  scheme 

uk  =  +  ^k+l  +  3»k+2 

corresponds  to  the  operator 
(43)  T  =  Y(E)  =  I  +  2E  +  3E2 
which  maps  O'  into  O'  ,  but  since 

Y(z)  =  1  +  2z  +  3z2 

has  the  roots  (-1  +  »J- 2  )/3  which  lie  in  the  unit  disc,  this 
map  is  not  always  pole-preserving. 

We  shall  now  state  and  prove  three  lemmas  which  collec¬ 
tively  serve  to  establish  the  above  theorem.  The  first 
focuses  on  the  banded  upper  triangular  structure  (35)  of  T, 
the  second  on  the  growth  condition  (37),  and  the  third  on 
the  requirement  that  Y(z)  be  zero  free  in  the  unit  disc. 

LEMMA  1 .  Let  T  be  a  linear  operator  which  maps  the 
space  of  transients  O'  into  itself.  The  following  are 
equival ent : 

(i)  c (E) r  y  =  O  whenever  y  €  7  and  c(z)  is  a 

polynomial  such  that  c (E)y  =0,  (i.e.,  if  y^ 

has  the  representation  (3)  for  some  choice  of  the 

a  ’s  and  z  ’s,  then  u,  =  T  y.  must  have  the  repre- 
v  v  k  •  k 
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1 


(U) 

(iii) 


sentation  (23)  for  some  choice  of  the  av's 
with  the  zv's  being  the  same.), 

T  has  the  representation  of  (32)— (35 )> 

T  commutes  with  the  shift  operator  E. 


Proof.  Assuming  (i)  we  see  that  since 

(E  -  .zl)  (1, z, z2, . . .  )  =  (0,0,0,...  ) 
we  must  also  have 

(E-zI )  r(l,z,z2,...  )  =  (0,0,0,...  ) 


for  each  choice  of  z  with  |z|  <  1  .  Using  the  general 
representation  (30)  for  the  linear  operator  T  we  find 


(E-zI)  r ( 1 ,  z ,  z  , . . .  ) 


CJO 


CO 


CO 


=  (E-zI )  (  ^  Y0^z  ,  ^2  yHz  •  SY2^Z  ^ 

£=0  1=0  £=0 


-(Yio+  £  (Y, 


£=1 


U  y0£-l)z 


co 

’  Y20+  ^  y2£~yi,£-l'>Z  ^ 

v  —  ± 


so  that  the  power  series 
oo 

Yk0  +  2  ^  Yk£-Yk-1,  £-l^z  ’  k=1’2"-- 

i=l 

must  vanish  for  every  choice  of  z  with  |z|  <  1  .  It  follows 
that 
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Y^q  -  0,  k  -  1,2,... 

Yk l  =  Yk-1 ,1-1  ’  k,t  =  » 

i.e.,  that  T  has  the  representation  (35)  or  equivalently 
(32)  so  that  (i)  implies  (ii). 

When  T  =  y(E)  has  the  representation  of  (32)-(33)  we 
clearly  have 

rEyk  =  ryk+l  =  Y0yk+l+Ylyk+2+Y2yk+3+.  .  .  =E  F  yk’ k=0’ 1  ’  •  •  • 

for  each  y  €  X  so  that  (ii)  implies  (iii).  Moreover,  if 
T  commutes  with  E,  then  T  also  commutes  with  the  polynomial 
c(E)  in  E  so  that 

c(E)  T  y  =  T  c(E)y  =  T  0  =  0 

whenever  c(E)  y  =  0,  i.e.,  (iii)  implies  (i  ).  (] 

LEMMA  2 .  Let  T  be  an  operator  having  the  representation 
of  (32)  -  (33).  In  order  that  T  map  CT  into  T  it  is  both 
necessary  and  sufficient  that  (37)  hold  (or  equivalently, 
that  (34)  be  analytic  on  the  disc  \z\  <  1.) 

Proof.  Assume  first  that  T  maps  7  into  T .  The  series 

2 

(34)  which  represents  the  first  component  of  T(l,z,z  ,...  ) 

2 

must  then  converge  whenever  (l,z,z  ,...  )  6  T  ,  i.e., 

whenever  |z|  <1  .  Thus  (34)  is  analytic  on  the  unit  disc 
so  that  the  Cauchy-Hadamar d  formula  (36)  must  produce  a 
radius  of  convergence  R  >  1,  and  (37)  holds. 
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Conversely,  assume  that  (37)  holds  and  that  y  €  O'  is 
given  so  that 

(44)  lira  sup  lyk|1//k  <  1  . 

Since  (44)  holds,  there  exist  constants  A>0,  0  <  &  <  1 
such  that 

(45)  |ykl  <  A  ak  ,  k=0,l,2f ...  , 

and  since  (37)  also  holds  there  exist  constants  B  >  0, 

0  <  3  <  a~^  such  that 

(46)  |Ykl  <  B  3k,  k=0, 1,2,...  . 

Using  (45  )  -  (46)  and  the  fact  that  0  <  a|5  <  1  ve  see  that 
the  series  (33)  which  is  used  for  the  k-th  component  of 
ry  is  majorized  by 


^ryk  I- I  YOyk+Ylyk  t-l  +  Y2yk+2+  ^ 


.  tjcO  .  A  k  nsl  .  .  k+1  n02  .  .  k+2 

<  Bp  Aa  +B3  Aa  +Bp  Aa  + 


=  ABak  /  (l-ag). 


and  is  thus  convergent,  k  -0,1,2,...  .  Thus  Ty  is  well 

defined,  and 

]  im  sup  |  ryR|1//k  <  lim  sup  |ABak/  ( 1  -ap  )  |  "*'^k-a  <  1 


so 


that  T  y  'J'  . 
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NOTE.  A  slight  extension  of  the  above  argument  shows 
that  the  operator  of  (32)  -  (33)  will  map 

Tr  =  |  (y0,y1,y2>  •  -  •  )T  :  lira  sup  |yk 1 1//k  <  R  ] 

into  itself  provided  that 

lim  sup  |  Yr  I  <  l/R> 

i.e.,  provided  \{z)  is  analytic  on  the  disc  \z\  <  R, 

0  <  R  <  cxd  .  The  case  R=1  is  covered  by  the  lemma. 

When  y(z)  is  analytic  on  the  unit  disc  and  |z|  <  1  ve 

have 

2  2 
y(E)  ( 1 ,  z  ,  z  ,...  )  =  Y(z)  (l,z,z  ,...  )  , 

and  ve  thereby  see  that  y(E)  will  annihilate  the  power 
sequence  associated  with  the  pole  z=zq  if  zq  is  a  root  of 
y(z).  More  generally,  suppose  that  y  has  the  representation 

co  , 

( 47 )  yk  =  ^  avzv  ,  k=0, 1,2,... 

v=l 

where 

co 

(48)  Uv|  <  cc 

v=l 

and  where 

(49)  |zv|  <  R  <  1  for  each  v=l,2,... 

^e  then  find 
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and  we  thereby  conclude  that  each  pole  zv  associated  with 
y  will  also  be  a  pole  of  a  =  y(E)y  provided  that  y(z)  has 
no  zeros  in  the  unit  disc. 


The  above  arguments  show  that  when  y(z)  is  analytic 
and  zero  free  on  the  unit  disc,  then  y(E)  is  a  pole  pre¬ 
serving  mapping  of  the  set  of  sequences  y  of  the  form  (47)- 
(49)  into  itself.  Although  such  transients  are  the  ones 
most  likely  to  be  met  in  practice,  for  the  sake  of  com¬ 
pleteness  we  shall  extend  the  argument  so  as  to  include 
the  somewhat  larger  class  0"  .  In  so  doing,  we  shall  find 
it  convenient  to  use  the  notation 

yU)  =  y0  +  yxz  +  y22  + 

for  the  generating  function  associated  with  a  given  sequence 
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y  =  iy0*yl»y2» • • •  ) »  with 

/  -lx  -1  -2 

yU  )  =  yQ  +  y2z  +  y2z  +  ... 

being  the  corresponding  z  -  transform.  Of  course,  when  y 
has  the  representation  (47)  -  (49)  we  find 

yu-1)=  y>-k  v  avzvk  = z  •  Eav/(z-zv) 

k=0  v=l  v=i 

so  that  z^,z2>...  are  the  poles  of  y(z  ^).  More  generally, 
we  shall  say  that  the  map  T  is  pole  preserving  on  X 
provided  that  for  every  choice  of  y  6  J  the  z  -  transform 
u(z_1)  of  u  =  Fy  has  exactly  the  same  poles  as  the  z-transform 
y(z-*)  of  y.  The  relationship  between  these  poles  is  made 
precise  in  the  following 

LEMMA  3 .  Let  y,  y  be  sequences  with 

(50)  lim  sup|y^|^^  =  R  <  1,  lim  sup|y^|^^  <  1 

and  let  u  =  y(E)  y,  ▼  =  y(E)Y  i  i.e., 

VVk,,lJk4l*Vk4!*'  •  •  ’  k=0>  1  ’ 2-  •  •  • 

WrtWVw*-'  k=o,i,2,...  . 

Then  the  generating  functions  v(z),  y(z)  are  analytic  for 
|z|  <  1,  the  z-transforms  y(z  ^),  u(z  ^)  are  analytic  for 
|z|  >  R,  and  the  identity 

u(z-^)  =  Y(z)*y(z-1)  -  v(z)  +  v(0) 


24 


holds  in  the  annulus  R  <  |z|  <  1  and  thus  in  the  common 
domain  of  analyticity  of  these  four  functions.  In  parti¬ 
cular,  when  y  is  zero  free  in  the  unit  disc,  the  poles  of 
u(z  *)  (which  must  lie  in  the  disc  |z|  <  R)  must  coincide 
with  the  poles  of  y(z-^). 


Proof.  By  using  the  Cauchy-Hadamard  formula  in 
conjunction  with  (50)  and  the  note  following  Lemma  2  we 
infer  that  y(z),  u(z)  are  analytic  for  |z|  <  l/R  and  that 
Y(z),  v(z)  are  analytic  for  |z|  <  1.  For  z  in  the  annulus 
R  <  |z|  <  1  we  then  have 

cc 


Y( z ) • y ( z  1)+v(0)=  Y,  E  Ykz*‘y£z  *+  E  Ykyk 


CO  CO 


-i 


o 

II 

o 

II 

X 

k=0 

-  E  wk~* 

+  E  \y^z 

-U-k) 

K>e 

e>k 

CO  CO 

cc 

oc 

-  E  ( 

v^)z  +  E  (  E  Ykyk+v)  z_ 

o 

II 

o 

II 

> 

v=0 

X 

ii 

o 

cc 

=  E  vvzV  + 

cc 

E  uv  2_v 

v-0 

v-0 

=  v  (  z  )  +  u  (  z 

.  0 

Taken  together,  these  three 

lemmas  give 

the  pr ev i ously 

stated  theorem  which  provides  a 

simple  chara 

c t  or i za t i on  for 

the  pole  preserving  mappings  we  would  like  to  use  in  a 
preprocessing  scheme.  The  discussion  of  section  2  can  now 


I 


» 

! 

i 


be  extended  at  once  to  the  case  where  the 
(38)  is  replaced  by  any  operator  T=  v(E) 
analytic  and  zero  free  on  the  unit  disc. 


> 


4 


'1 


tail  sum  operator 
for  which  y(z)  is 


26 


4.  NOISE  SUPPRESSION 


Let  a  ^  0  and  z  be  given  with  |z|  <  1,  and  let 

(51)  vk  =  azk  +  ek,  k=0,l,2,... 

where  eo>ep>£2’***  are  independent  random  variables  with 
common  mean 

<  ek  >  =  0,  k=0,l,2,... 
and  common  variance 

<  ek  >  =  a2,  k=0,l ,2, ...  . 

Let 

F  =  y(E)  =  y0+Y1E-tY2E2+  ... 
be  a  preprocessing  operator  for  which 

(52)  Y2  =  1 Y0 1 2  +  |Yj  I2  +  |y2|2  +  ... 

is  finite,  and  let  ▼  =  IV  so  that 

(53)  vk  =  ay(z)zk  +  6R,  k=0,l,2,... 
where 

(54)  6k=Y(E)€k=Y0£ll+Y1ek+1  +  Y2ek+2  +. . . ,  k=0,l,2, 
is  a  random  variable  with  mean 

CC  CC 

(55)  <6k>  .  <  E  =  E  vw°>  k=on.2,... 

i= 0 
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and  variance 


2  o°  _  co 

(56)  <|6ki  >=<  E  v4Vt+tW=  EVAn'kV 

t  >  1-1=0  l ,  (j=0 

oo 

=  £  °2=Y2a2,  k=0, 1,2,... 

£=0 

The  processed  noise  6 q,  6^,  62>'**  i s  correlated  with 

oo 

(57)  <6k6k+P>=  £  ',Vw  ek+P+M 

l,  U=0 

■<VoVl  +  Yp+2Y2+*  •  *  »  k  ,  p-0, 1,2,... 

Ve  would  like  to  develop  some  quantitative  measure  of 
the  tendency  of  the  preprocessing  scheme  to  suppress  the 

i  k  i  2 

noise.  Ve  shall  use  |az  |  as  a  measure  of  the  signal 

•  7  2 

present  in  the  k-th  component  of  w  with  =  o"~  beirtg 

a  corresponding  measure  of  the  noise.  Analogously, 

|ay(z)z^|^  gives  a  measure  of  the  signal  in  the  k-th 

component  of  v  =  Tw  with  <|dk|2>  =  y2o2  being  the 

corresponding  measure  of  the  noise.  The  preprocessing 

scheme  r  thus  improves  the  k-th  component  s i  gna 1 -to-noi se 

ratio  by  the  factor 

-  iyu)i2/y2 

|az  |  /o 

Since  this  ratio  is  independent  of  k,  ve  see  that 
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2 


(58)  g(z)  =  |y(z)|2/y2  =  |y0+y1z+y2z2+  ...  1 

|y0I2+Iy1I2+|y2|2+..^ 

provides  us  with  a  quantitative  measure  of  the  tendancy  of 
r  to  enhance  the  signal-to-noi se  ratio  for  a  power 
sequence  based  on  the  pole  z  .  When  g(z)  >  1,  the  noise 
is  suppressed,  and  we  refer  to  the  set 

(59)  =  \z  €  C  :  |z|  <  1  and  g(z)  >  l| 

as  the  region  of  pole  amplification  associated  with 
r  =  Y(E). 

Suppose  now  that 
n  , 

(60)  w^  =  Y  avzv  +  ek’  k=0,l,2,... 

v=l 

where  a^O  and  |  zv  |  <  1  for  v=l,  .  .  .  ,  n,  where  z^,  .  .  .  ,  zq 

are  distinct,  and  where  6^,0^, e 2,...  are  again  independent 

2 

random  variables  with  zero  mean  and  common  variance  o 
If  we  compute 

(61)  vk  =  y(E)wk  =  Y  avY^zv^zvk+  k=0»1>2»--- 

v=l 

(with  ^o,^l,^2’**'  aSa^n  given  by  (54)  )  we  will  expect  to 

enhance  the  pole  zv  relative  to  the  noise  provided  that 

zv  lies  within  the  region  of  amplification  Cl  with  the 

degree  of  enhancement  depending  on  the  size  of  g(zv).  In 

practice,  the  poles  z,....,  z  are  unknown,  but  we  often 

I  n 
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have  some  a  priori  knowledge  of  the  portion  of  the  unit 
disc  in  which  they  are  most  likely  to  be  found.  In  such 
situations  we  might  reasonably  select  a  preprocessing 
scheme  which  is  specifically  designed  for  signals  of  the 
type  we  expect  to  process,  i.e.,  we  choose  r  =  y(E)  so  as 
to  make  g(z)  as  large  as  possible  in  the  region  of  the 
unit  circle  where  we  expect  to  find  z^,...,  zn*  An  upper 
bound  on  the  possible  size  of  g(z)  is  provided  by 

THEOREM  2 .  Let  the  operator  r  =  y(E)  satisfy  (32)  - 
(34),  (37)  v’ith  (52)  being  finite  and  let  g(z)  be  the 
corresponding  signal-to-noi se  amplification  factor  (58). 
Then 

(62)  g(zQ)  <  l/(l-|z0|2),  I  z0 1  <  1 

with  equality  holding  if  and  only  if 

(63)  Y(z)  =  a/ ( 1 -ZqZ ) ,  |zj  <  1 

where  a  is  a  nonzero  constant  so  that  y(z)  has  a  simple 
pole  at  the  point  of  inversion  z  =  1/z"q  of  Zq  relative  to 
the  unit  circle  and 

(64)  Tyk  =  a(yk+z0.yk+1+z02.yk+2+  ...  ),  k=0,l,2,...  . 

Proof.  Using  Cauchy's  inequality  we  find 
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(65)  |y(z) I2  =  j  £  Yk  zk j2 

k=0 

<  E  l\l2  E  Ukl2 

k=0  k=0 

=  Y2/ ( 1 - | z | 2 ) 

vith  equality  holding  if  and  only  if  y  is  a  scalar  multiple 
of  the  pover  sequence  (l,z,z  ,...  ).  In  conjunction  vith 

(58)  this  gives  (62)  -  (64).  § 

The  bound  (62)  provides  a  natural  limit  to  the  signal- 

to-noise  amplification  ve  can  achieve  tnrough  the  use  of 

a  preprocessing  scheme.  'when  the  sampling  interval  T  is 

chosen  so  large  that  the  pole  zv  lies  deep  vithin  the  unit 

o  3 

circle,  i.e.,  zv~  0 >  the  successive  povers  zv  ,  z  v 
decay  so  rapidly  that  substantial  noise  suppression  is 
impossible  and  ve  find 

g(zv)  <  1/(1- |  2v  | 2 )  =;  1  . 

On  the  other  hand,  if  T  is  so  small  that  zsi  lies  near  the 
rim  of  the  unit  circle,  i.e.,  |zv|  <  1  but  |zv|  z:  1,  then 
the  maximum  of 

g(zv)  =  ]/(l-|zv|2)  ^  1  /  [  2  ( 1  —  |  z.y  1  )] 

VI 11  be  large  and  a  suitable  preprocessing  scheme  can 
achieve  a  substantial  suppression  of  the  noise,  cf.  Fig.l. 
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Suppose  now  that  we  have  some  a  priori  knowledge  tnat 
the  poles  z^,...,  z which  we  seek  all  lie  in  some  subset 

CL  of  the  unit  disc,  e.g.,  if  the  sampl ing  interval  T  is 

small  we  expect  to  find  z^,...,  z^  in  a  tight  cluster  near 
z=l  so  we  might  take  Q.  to  be  the  lens  shaped  set 

(66)  CL  =  { z  €  C  :  |z|  <  1  and  |  z  —  1  |  <  pi 

obtained  by  intersecting  the  unit  disc  with  the  disc  of 
radius  p  which  is  centered  at  z=l.  Once  Ol  is  chosen,  we 
would  like  to  determine  a  corresponding  preprocessing 
operator  r  -  y(E)  which  is  optimal  in  the  sense  that  the 
minimum  value  taken  by  g(z)  as  z  ranges  over  CL  (i.e.,  the 
smallest  s i gnal - to-no i s e  amplification  at  any  pole  location 
we  might  possibly  encounter)  is  as  large  as  possible.  Vnen 
CL  ~\zq\>  ^he  optimal  operator  r  is  given  by  (63)  -  (64). 

In  more  realistic  situations  such  as  (66),  tnere  is  no 
known  procedure  for  constructing  an  optimal  r.  Nevertheless, 
we  can  develop  good  if  not  best  preprocessors  for  certain 
natural  choices  of  CL  and  we  shall  now  proceed  to  show  how 
this  is  done. 

Suppose  first  that  Cl  is  given  by  (66)  with  p  being  small. 

V'nen  y(z)  is  given  by  (39)  the  contour  lines  of  |y(z)|  are 

2  _  2 

circles  centered  at  z=l  with  |y(z)|  >  P_  when  0  <  | z  —  1 1 <  p. 

Ve  might  thus  expect  the  tail  sum  operator  S  =  y(E)  of  (38) 
to  be  an  ideal  choice.  Unfortunately,  the  corresponding 
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so  we 


2 

Y  of  (52)  is  not  finite  when  y(z)  is  given  by  (39), 
are  forced  to  consider  related  approximations. 

One  very  natural  approximation  is  obtained  by  simply 
truncating  the  infinite  series  (39)  so  as  to  obtain 
(42)  with  N  being  a  positive  integer.  The  corresponding 
preprocessing  scheme  u  =  y(E)y  is  then  the  local  averaging 
procedure  of  (40),  and  the  si gna 1-to-noi se  enhancement 
factor  is  given  by 

(67)  g(z)  =  1 1  +  z+  ...  +z^-1|2/N  =  |  l-z^  |2/(  1 1-z  |2N)  . 

By  using  Cauchy's  inequality  we  see  that  this  operator  is 
optimal  in  the  sense  that  it  maximizes 

g(l)  =  1  W+  •"  +VN-1Z'  I 

Iy0I2+I>i  l2+  •  •  •  +Iyn_1  I2 

as  y(z)  ranges  over  the  set  of  all  (nonzero)  polynomials 
having  degree  less  than  N.  If  we  let 

Gk  =  exp ( i  •  2nk/N),  k=l,2,...,  N-l,  (i2=  -l) 

denote  the  N-th  roots  of  unity  other  lhan  z=l ,  we  may  write 
(67)  in  the  form 

g(z)  =  |z-c1l2  ...  |z-Cn_1I2/n 

and  tnereby  conclude  that  the  contour  lines  of  g(z)  coincide 
with  the  equipotenti a  1 s  of  the  two  dimensional  field  wnich 
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results  when  unit  charges  are  placed  at  z  =  ^ i ’ ^2 5  *  *  *  ^  j 
in  the  plane.  This  helps  us  visualize  the  region  of 
amplification  when  N  is  small,  cf.  Fig.  2.  For  large  N 
we  neglect  z^  in  (67)  to  see  that  the  region  of  amplifi¬ 
cation  is  essentially  the  lens  shaped  set  (66)  which  results 

_l/2 

when  we  take  p  =  N  and  that  g(z)  increases  from  approx¬ 

imately  1  to  its  maximum  value  N  as  the  distance  from  z  to 

_  1  /  2 

1  decreases  from  N  '  to  0.  In  this  way  we  see  that  a 
large  signal-to-noise  amplification  is  possible  only  when 
N  is  large  in  which  case  the  corresponding  region  of  pole 
amplification  is  small.  Thus  N  should  be  chosen  in 
conjunction  with  the  sampling  interval  T,  i.e.,  when  T  is 
small  the  zv's  are  tightly  clustered  near  z=l  and  a  large 
value  of  N  is  appr opr  rate,  cf.  Fig.  1  and  Fig.  3 

A  second  natural  approximation  to  (39)  is  obtained  b}1- 
slightly  shifting  the  pole  z=l  of  y  to  the  nearby  point 
z  =  R  >  1  so  as  to  make 

(68)  y(z)  =  1  +  z/R  +  z2/R2  +  ...  =  (l-z/'R)"1 
In  this  case  u  =  y(E)y  is  given  by 

(69)  uk  =  yk  +  yk+1/R  +  >'k+2/R2  +  •  •  ■ 

and 

(70)  g(z)  =  (R2-l)/|z-R|2  . 

The  contour  lines  of  g(z)  are  circles  centered  at  z=R 


34 


with  g ( 2  )  =  K  on  tJje  circle  vnere 
|z-R|  =  [ (R2-l )/K]l/2  . 

In  particular,  the  maximum  signal-to-noise  enhancement  is 
g(l)  =  (R+l )/(R-l )  s  2/ (R-l ) 

and  the  region  of  amplification  is  that  portion  of  the  unit 
disc  jz|  <  1  which  lies  within  the  circle 

|z-R|  =  (r2-i)1//2  . 

Again  we  observe  a  fundamental  trade  off  between  the  size 
of  the  region  of  amplification  and  the  maximum  possible  g, 
c  f .  Fig.  4 . 
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CONSTRUCTING  NOISE  SUPPRESSING  PREPROCESSING  OPERATORS 


P  . 


In  many  cases  of  practical  importance  we  encounter  a 
system  which  exhibits  various  modes  of  oscillation  with 
frequencies  which  are  more  or  less  integral  multiples  of 
some  fundamental  frequency  and  with  all  of  these  modes 
being  similarly  damped  (cf.  Appendix  l),  i.e.,  the  s^'s 
of  (1)  are  more  or  less  uniformly  spaced  along  some  line 
Re  s  =  -a  in  the  left  half  plane  and  the  corresponding  zv '  s 
of  (4)  are  more  or  less  uniformly  spaced  around  the  circle 
|z|  =  exp(-aT)  <  1  .  To  aid  us  in  the  numerical  extraction 
of  the  first  few  zv's  we  would  like  to  use  a  preprocessing 
scheme  which  is  optimal  on  some  corresponding  sector 

(71)  Q.  =  |z=rei6  :  p  <  r  <  1,  -Q  <e<0l 

of  some  annulus  of  the  unit  circle  with  tne  parameters 
p  and  0  (0  <  P  <  1  and  0  <  0  <  n)  being  determined  by  the 
damping  coefficient  a,  the  fundamental  frequency,  the 
sampling  interval,  and  the  number  of  poles  we  are  trying 
to  find.  Tne  preprocessing  schemes  of  (40)  and  (69)  (which 
are  designed  for  the  region  (66)  )  tend  to  enhance  the 
lowest  modes  but  damp  out  any  others.  We  must  thus  devise 
a  preprocessing  scheme  r  =  y(E)  for  a  region  of  the  form(7l). 
Ve  now  describe  one  way  in  wnich  this  can  be  done. 

We  first  observe  that  the  function 
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(72)  Y*U)  = 


z^-2Rz  c  os  6  +R^ 


_ R 

(z-He^)(z-Re”^) 


is  analytic  inside  the  circle  |z|  =  R  on  vhich  its  tvo 
conjugate  poles  z  =  R  exp(  +  i0)  are  found.  For  |z|  <  R 
and  0  <  9  <  rt  we  may  write 


Y*(z)  = 


i  0 


-i0 


2i  sin0  (  1-  (  z/R )  e 


10 


l-(z/R)e  16  ) 


1 

2i  sin  0 


£  (el6z/R)k-e-ie  £(e-i6z/R)k', 
k=0  k=0 


and  after  a  bit  of  simplification  this  yields  the  series 
r epresenta ti on 

2 

iTii  \  i  ?m.3  z  sin  z 

Kti)  Y*(z;  -  1  +  - -  •  —  1-  - -  •  — p  +  ... 

sin0  R  sin  6  h 


for  {72). 

To  obtain  our  desired  preprocessing  operator  we  shall 

now  select  a  number  p  such  R  ,  6.  pairs  with  R  >  1  and 

J  3  3 

0  <  0.  <  tt ,  and  then  form  the  product  of  the  corresponding 
3 

terms  (72),  i.e.,  ve  take 


R 


(74)  Y(z)  =  ff  -j 


.  z  -2R . z  cos  0  .  +R  .  “ 
3-1  3  33 


2  =  Y^+YiZ+Y0z 


'0  '1 


From  (73)  -  (M)  ve  see  that  y(z)  has  the  alternative 
representation 


37 


,JL-  f  sin2 0  . 

(75)  y(z)  =  TT  I1  +  - 1  ~ 

,  l  sin  9  .  R 

0  =  1  0 


sin36  •  z‘ 


sin  9.  R.  sin  8.  R  . ‘ 

JO  J  J 


...  |  > 


so  that  if  we  define 


Y(1^(z)  =  1  + 


sin20,  z  sin39,  z‘ 


sin  9^  R^  sin  9^  R^‘ 


+  .  . 


and  successively  compute 

r«>U)  .  v<‘-1)(z).{  i+  ^ 

sin  e£  R£  sin  e£  Re  j 


for  1= 2,3,...,  p  the  resulting  y'P'(z)  will  coincide  with 
(75)  and  thus  (74).  Ve  may  thus  numerically  obtain  the 
coefficients  Y q,Y j,Y 2>  • • «  °f  (74)  by  using  the  recursion 
scheme 


, ( 1 ) _  sin(k+l )  9,.  _  1 

k  ~  sin  8,  "pi 


,  k=0,l,2. 


.!“■  e  .  y  ■  >."..,1 . 

j=0  sin  ne  1  =  2,3,  ...  ,j> 


with 


(77)  Yk  =  YkP^,  k=0, 1,2, 


Once  Yq>Y^^Y2>...  are  known,  we  use  the  corresponding 


operator 
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r  =  Y0  +  VjE  4  v2E  4  ... 


for  our  preprocessing  scheme,  i.e.,  we  take 
uk  =  YOyk  +  Ylyk4l  +  Y2yk42  +  ’ * ‘ 

Prom  (74)  we  see  that  y(z)  will  be  analytic  on  the  disc 
|  z  |  <  R  when 

R  =  min  {R^,...,  R  1  >  1  . 

By  using  the  Cauchy-Hadamard  formula  (36)  we  then  see  that 
there  are  constants  C>0,0<q<l  such  that  the 
coefficients  Yq, Y2> • • •  ve  compute  in  (79)-(80)  will 
satisfy  the  bound 

I  Yi^  |  <  C  •  q  ,  k=0, 1,2,..,  . 


Ve  thereby  infer  that  the  sum 

E  l\|2  <  E  |Cqk|2  =  C2/(l-q2) 


k=0 


k=0 


is  finite  so  that  the  corresponding  g(z)  of  (58)  is  well 

defined.  Our  goal  is  to  determine  the  pole  location 

parameters  R.,0.,  j=l,...,  p  which  result  in  a  g(z)  which 
3  3 

is  uniformly  large  on  all  of  &  in  the  sense  that 


min 


|  g(z)  :  z  €  &  | 


is  as  large  as  possible.  Since  g(z)  is  the  modulus  of  an 
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analytic  function  which  is  zero  free  on  (X  ,  this  is 
equivalent  to  maximizing  the  minimum  value  taken  by  g(z) 
as  z  ranges  over  the  boundary  of  Q.  . 

Using  a  digital  computer  in  an  interactive  mode 
(where  it  is  possible  to  perturb  the  pole  location  parameters 
R  ,  6.  and  immediately  see  the  effect  on  the  contour  lines 
of  g(z)  we  have  been  able  to  design  good  if  not  best 
preprocessing  operators  T  =  y(E)  for  the  region  (71),  cf. 
Fig.  5.  More  work  is  needed  in  order  to  systematize  this 
procedure,  however. 


40 


6.  ALTERNATIVE  METHODS  FOR  COMPUTING  THE  NULL  VECTOR 


I 


At  the  heart  of  transient  analysis  lies  the  problem 
of  extracting  a  suitable  null  vector  c  =  (cq,c^,...,  cfl) 
for  the  data  matrix  ^  of  ( 1 1 ) .  As  we  have  seen,  c  can 
be  obtained  alternatively  as  a  null  vector  of  certain 
related  matrices  such  as  (25),  (26)  or,  more  generally, 
as  a  null  vector  of  the  corresponding  matrices  which  result 
when  the  tail  sum  operator  S  is  replaced  by  some  other 
preprocessing  operator  r  chosen  for  its  noise  suppressing 
properties.  Ve  shall  now  describe  a  general  procedure  of 
this  type  in  the  hopes  of  obtaining  a  more  robust  scheme 
for  computing  c  and  thereby  for  computing  z.,  .  .  .  ,  . 

Indeed,  suppose  that  we  are  given  the  transient 
sequence  y  =  (Yq)  y  ^  >  yp »  *  *  •  )  an^  f*16  preprocessing 

operators  rv  =  Y^(E),  i=0,l,...,  m.  Ve  shall  numerically 
generate  the  sequences  P^y,  i=0,l,...,  m  and  then 
assemble  the  matrix 

fVo  1  0yl  '  •  •  ^0yn 

(78)  Jtj  =  fVo  riy2  •  •  •  riyn 

^"my0  ^myl  ’  '  ’  ^myn 


Using  a  singular  value  decomposition  of  or  an  eigenvalue 

T 

decomposition  of  ^  as  described  in  section  1  we  shall 
compute  a  null  vector  c  =  (cq,c^,...,  cr)  for  and  then 
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take  z. , .  .  .  ,  z  to  be  the  roots  which  result  from  the 

corresponding  factorization  (6).  Of  course,  tnis  scheme 

will  be  successful  only  if  the  operators  are  chosen  in 

such  a  manner  that^ft  has  the  same  null  space  as  the  data 

matrix  y  .  This  imposes  certain  constraints  on  the 

operators  T,-.,  P,  ,  P  which  we  shall  now  explore  in 

Ulm 

some  detail. 

Let  Yq(z),  Y^(z),...,  Ym(z)  be  analytic  on  the  unit 
disc  J z |  <  1  We  shall  say  that  these  functions  satisfy 
the  Haar  condition  provided  that  no  linear  combination 

(79)  Y„U)  =  b0Y0U)  +  bjY,(2)  +  ...  +  bm  Y.IO 

of  these  functions  lias  more  than  m  zeros  in  |z|  <  1  units 
b.,,  bn  ,  b  all  vanish.  After  illustrating  this 

concept  with  several  examples  we  shall  show  that  this  is 
exactly  the  property  which  is  needed  in  order  to  insure 
that^  and  y  share  the  same  null  space. 

Example  1.  The  functions 

Y1  ( z )  =  z1 ,  i=0, 1  ,  .  .  .  ,  m 

satisfy  the  Haar  condition  on  the  unit  disc  since  (79) 
reduces  to  the  polynomial 

Y* ( z )  =  b0  +  blZ  +  ...  +  bmzm 

which  has  at  most  m  zeros  unless  b„  =  b.  =  . . .  -  b  =0. 

0  1  m 

In  this  case 
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r  y  = 


y3  = 


yi+j 


so  that^j  is  the  matrix  which  results  when  all  but  the 
first  m+1  rows  of  are  deleted. 

Example  2.  When  Y(z)  is  analytic  and  zero  free  on 
the  unit  disc  the  functions 


(80)  Y^z)  =  zi  Y  (  z  )  ,  i=0, 1 , .  .  .  ,  m 

satisfy  the  Haar  condition  there  since  (79)  reduces  to  the 
product 

Y,U)  =  Y(z)  •  [bQ  +  V  +  ...  +  b^”] 

which  has  at  most  m  zeros  unless  b„=b-.  =  ...  =  b  =0.  In 

0  1  m 

this  case 

TVj  =  EiY(E)y  .  =  Y(E)yi+j 


so  tnat^j  is  obtained  by  generalizing  (25)  to  the  case 
where  is  replaced  by  y(E). 

Example  3.  Again  let  y(z)  be  analytic  on  the  unit 
disc  and  assume  xhat  y(z)  is  also  one-to-one,  i.e.,  that 
Y(z)  =  y(z')  only  when  z  =  z'  .  If 

( 8 1 )  Yi(z)  =  y(z)A  ,  i=0, 1 , . . . ,  m 

then  (79)  reduces  to  the  polynomial 

Y*(z)  =  bQ  +bj y( z )  +  b2Y(z)2  4  ...  4  bmY(z)m 
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of  degree  m  in  v(z).  Since  y(z)  is  one-to-one,  y*  ( z )  will 

have  at  most  m  zeros  unless  b„  =  b.  =  . . .  =  b  -  0,  and 

0  1  m 

so  the  Haar  condition  is  satisfied  once  again.  In  this  case 

riyj-v(E)iyj 

so  that  is  obtained  by  generalizing  (26)  to  the  case 
where  S  is  replaced  by  y(E). 

Example  4.  Let  £q>  £^,...,  £m  be  distinct  complex 
numbers  with  I I  <  1  for  each  i,  and  let 

(82)  y.(z)  =  1/(1-  Q1  z)  ,  i=0,l,...,  m  . 

Ve  may  express  the  corresponding  expression  (79)  in  the  form 
Y*(z)  =  b0/(l-C0z)  +  ...  +  bm/(l-Cmz) 

=  d0  +  diz  +  +  <Vm 

(l-£0z)(l-C1z)...(l-Cmz) 

with  d^,  d^,...,  d^  depending  linearly  on  b^bj,...,  bm> 
Clearly  Y*.(z)  has  at  most  m  zeros  unless  dn=d1=  ...  =  d  =  0 
in  which  case  bQ=b^=  ...  =bm=0  also,  so  the  y^(z)'s  satisfy 
the  Haar  condition  on  the  unit  disc.  Using  Theorem  2  we 
see  that  the  operator 

P.  =  y.(E)  =  I  +  C-E  +  C^E2  +  .  .  . 

provides  the  maximum  possible  improvement  in  the  signal- 
to-noise  ratio  for  a  pole  at  z=£^,  and  thus  such  a  pole  will 
be  strongly  represented  in  the  i-th  row  of ,  i=0, !,...»  m. 
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Using  thiSy^J  to  determine  the  null  vector  c  may  be  preferable 
to  using  (25),  (26),  or  the  generalization  of  Examples  2,  3 
above  when  the  system  poles  are  veil  separated  and  it  is 
difficult  or  impossible  to  find  a  single  operator  P  which 
has  a  region  of  amplification  containing  all  of  these  poles. 

Example  5.  Again  let  £q>  £^,...,  £m  be  distinct 
complex  numbers  vith|£^|  <  1  for  each  i  and  let 

Y-(z)  =  TT  “ — ^  ,  i=0.  It...,  m 
J*  Ci‘Cl 

be  the  Lagrange  interpolating  polynomials  of  degree  m  which 
are  based  on  the  m+1  points  Cq>  £^,...,£ra  so  that 

Yi(Cj)  =(l  if  i= j  =  0,1,...,  m 

to  if  irj  . 


Since  YQ(z),  Y^(z),  ...,  Ym(z)  are  linearly  independent 

polynomials  of  degree  m,  the  Haar  condition  is  satisfied. 

The  application  of  the  operator  P^  =  Y^(E)  to  a  given 

sequence  y  =  (y^y^,...  )  can  be  effected  by  successively 

applying  the  m  operators  E  -  C ^ I ,  j^i  ,  and  then 

suitably  scaling  the  resulting  sequence. 

We  note  that  in  the  special  case  where  y  has  the 

representation  (3),  where  m=n-l,  and  where  we  take 

C_=z, ,  C,=z~,...,  C  -z  ,  we  find  that  the  elements  of 
0  11  2  m  n 

Jb  are  given  by 
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fVj  =  E 

v=l 


avvi<zv>2vJ  =  *i+1 


Thus  ^  is  a  row  scaled  version  of  the  Vandermonde  matrix 

based  on  z  ,  z 

l  n 

THEOREM  3 .  Let  ("V  =  Y^(E)  where  Y^(z)  is  analytic  for 
| z |  <  1,  i=0,l,...,  m  and  assume  that  Yq(z),  ( z Ym ( z 

satisfy  the  Haar  condition  on  tne  unit  disc.  Let 
y  =  (y0,y1,y2# • • •  )  be  a  transient  sequence  having  the 
representation  (3)  with  av^0  and  jz^J  <  1  for  each  v=l,...,  n 
with  z^,...,  zn  being  distinct,  and  with  n  <  m+1.  Then 
c  =  (cQ,c  , . . . ,  cq)  is  a  null  vector  of  the  matrix  >u  of 
(78)  if  and  only  if  c  is  also  a  null  vector  of  the  data 
matrix  of  ( 11 )  . 

Proof.  Using  (3)  we  see  that  for  any  choice  of 

cn,  c ,  c  we  have 
U  1  n 


(83)  c0^0+c1^il  +  ...+cn^in=(c0+c1E-f...+cnEn)Yi(E)  £avz, 

v=l 


k=0 


n 

=  £av(  co+cizv+'  •  *+cnzvn^Vi  ^zv^* 

V=1 


Now  if  c  is  a  null  vector  of  ,  the  factorization  (6)  must 
hold  so 

cq  ^3  j  o+C  1”^  i  1 +  *  ’  *  +Cn  ^  i  n  =  m 


and  thus  c  is  also  a  null  vector  o  fh  .  Conversely,  if  c 


is  a  null  vector  oi'Jti  the  left  side  of  (83)  vanishes  for 
each  i=0,l,...,  m  and  therefore 


(84) 


n 

£a 

v=l 


v(cO+clzv+‘ 


+cnZv 


n)[Vo( 


zv)+blYl(zv)+---+bmYm(zv)] 


for  every  choice  of  bQ,  b^...,  b^.  Since  yQ(  z ) ,  y  (  z ) , .  .  . , 
Ym(z)  satisfy  the  Haar  condition,  we  can  choose  the  b^ ' s 
so  as  to  make  (79)  interpolate  any  given  m+1  points,  and 
in  particular  (since  m  >  n-l)  we  can  choose  these  coefficients 
so  as  to  make 


(85)  b0Y0(zv)+blVzv)  +  --'+bmYm(zv> 


=  av(cO+ClZv+* • -+CnZvn)’  v=1"*-’  n 


Upon  substituting  (85)  into  (84)  ve  find 

n  2 

£  K(c0+ClZv+-  •  ‘+CnZvn)  1  =  0 

v=l 

so  that 


cny . +  c , y . , , +. . . +c  v . 

Cr  l  l*7  i  +  I  n"  l+n 


r  £av<c0+QlZv+* 
V=1 


•  •+CnZvn)zv1=0’1=0--1> 


i.e.,  c  is  a  null  vector  of  . 


The  ideas  leading  to  the  construction  of  (28)  lead  to 
the  somewhat  more  general  result  which  contains  Theorem  3 
as  a  special  case. 


THEOREM  4.  Let  r,i=  V^E)  where  Y^z)  is  analytic  for 

jz[  <  1  ,  i=0,l,...,  m  and  assume  that  Yq(z),  ( z ) ,  ...» 

Y  (z)  satisfy  the  Haar  condition  on  the  unit  disc.  Let 
m 

A-  =  6(E)  where  6.(z)  is  analytic  for  |z|  <  1,  j=0,l,...n 
J  J  J 

and  assume  that  6q(z),  6^(z),  ...,  6n(z)  also  satisfy  the 
Haar  condition  on  the  unit  disc.  Let  y  be  a  transient 
sequence  having  the  representation  (3)  with  a^/  0  and 
|z^|  <  1  for  v=l,...,n  ,  with  z^,...,  z^  being  distinct, 
and  with  n  <  m+1.  Then  d=(d^, d ^ , .  .  .  ,  d^)  is  a  null 
vector  of  the  (m+l)  x  (n+l)  matrix 


roAoyo  ^0Aly0  '  *  '  Wo 
^iAoyo  r* lAlyO  •  •  •  r iAny0 


(86)  7 n= 


^mVo  rmAlyO  *  *  '  PmV( 


if  and  only  if  the  corresponding  characteristic  function 

(87)  q ( z )  =  dQ  6q(z)  +  d  1 ( z  )  +  ...  +  d^iz) 

vanishes  at  z. , . . .  ,  z  and  is  elsewhere  nonzero  for  z  <  1, 
I  n 

Proof .  Using  che  representation  (3)  we  see  that 


(88)  bmd=  Eav[b0Y0(zv)  +  b1Y1(zv)  +  .  .  .+bmYm(zv)] 
v=l 

td060(zv)+dl6l(zv)  +  ---+dn6nUv^ 
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T  T 

for  every  choice  of  b  =(bQ,b^,...,  bm )  and  d  =  (dQ, , . . . , ) . 

If  d  is  chosen  so  that  (87)  has  z,  , . .  .  ,  z  as  roots,  then 

1  n 

the  left  side  of  (88)  must  vanish  for  every  choice  of  b 

(and  in  particular  for  the  choice  b  =  7nd)>  and  we  conclude 

that  d  is  a  null  vector  of  IT)  • 

Conversely,  suppose  that  d^O  is  a  null  vector  of  7)1 

so  that  the  right  hand  side  of  (88)  vanishes  for  each 

choice  of  b.  Since  yAz),  Y-,(z),...,  Y  (z)  satisfy  the 

u  l  m 

Haar  condition  and  n  <  m+1  it  is  possible  to  choose  b  so 
that 

l,0T0(2>)+blYl(2»)  +  ---+bmYn.(2») 

=«v[d060(2»>+dl6l(2v)  +  '--+dn6r,(2v)] 
when  v=l,...,  n.  It  then  follows  that 

£lav[W2v>+dlV2v>  +  ---+d„6„<2v>]  I2  -  0 

V=1 

and  thus  (87)  must  vanish  at  z, , . . . ,  z  .  Moreover,  since 

i  n 

6q(z),  6(z),...,  6n(z)  also  satisfy  the  Haar  condition,  the 
characteristic  function  (87)  has  no  other  zeros  in  the 
unit  disc.  [] 

Example  6.  Let 

Yi(z)  =  z1,  i =0 , 1 , . . . ,  m 

and  let 
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1 


6j(z)  =  Y(z)J,  j  =0,1,...,  n 

where  y(z)  is  analytic  and  one-to-one  for  |z|  <  1  .  As 

noted  in  Examples  1,3  above  the  y.(z)'s  and  6.(z)'s  satisfy 

^  J 

the  Haar  condition  on  the  unit  disc.  If  we  set  P  =  y(E) 
we  may  write 


Px  =  Y^E)  =  E1,  i=0, 1 ,  .  .  .  , 


m 


A j  =  6^(E)  =  P3,  j=0, 1 , . . . ,  n 


so  that 


(89)  TD  = 


0 


ry 


o 


,—1  n 

.  P  y 


o 


n 


ri  ry i  ...  r  yi 


ym  pym  ■  •  •  p"y„ 


Ve  then  conclude  that  d  =  (d^,  d^,...,  d^)  is  a  null  vector 


of  TT]  if  and  only  if  y(z^),...,  Y(zn)  are  roots  of 


the 


polynomial  d^  +  d^w  +...+  dnw  ,  and  since  Y  is  one-to-one 
the  z^'s  are  uniquely  determined  from  d.  Ve  observe  that 
(28)  is  the  specialization  of  (89)  to  the  case  where 
Y ( z )  =  l/(l-z)  as  used  in  Jain's  analysis. 

Example  7.  Let 


Yi ( z)  =  z1 ,  i=0, 1, . . . ,  m 


50 


and  let  6q(z),  6^(z),...,  6n(z)  be  the  Chebyshev  polynomials 
which  may  be  defined  and  computed  by  means  of  the  three 
term  recursion 

60(z)  =  1 
(90)  §1(z)  =  z 

6.(z)  =  2z6  •  , ( z )  -  6  2(z),  j=2,3,...,  n  . 

J  J  J 

Since  6(z)  is  a  polynomial  of  exact  degree  j,  the  Haar 
J 

condition  is  clearly  satisfied.  By  using  (90)  we  see  that 
it  is  possible  to  numerically  generate  the  sequences 

A.y  =  6.(E)y  ,  j=0,l,...,  n 
J  J 

when  y  is  given  by  successively  computing 

Vi  =  yi’  . 

Vi  =  yi+i»  i=0»1»2»-“ 

Ajyi  =  2Aj-lyi  +  l  "  A j -2y i *  1=0> 1 >  2 , . . . ,  j=2»3,...,  n. 

Vhen  y  has  the  representation  (3)  and  m  >  n-1,  z^,...,zn 

will  be  the  roots  of  the  characteristic  polynomial  (87) 

(which  is  now  parametrized  using  the  Chebyshev  polynomials) 

T 

if  and  only  if  d  =  (dQ,d^,...,  dR)  is  a  null  vector  of  the 
matrix 

A0y0  Aly0  •  •  •  Any0 

Vl  Alyl  •  •  •  anyx 

•  •  • 

A_y  A,v  Ay 

m  1  “  m  n^  m 
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7.  THE  PROBLEM  OF  EXTRANEOUS  ROOTS 


If  in  the  process  of  numerically  analyzing  a  given 
transient  y  =  (yQ/yj#y2**-  )  of  the  form  (3)  we  overestimate 
the  system  order  n,  i.e.,  we  carry  out  the  calculations  with 
n  replaced  by  some  larger  integer  m=n+r>  we  obtain  a 
characteristic  polynomial 

(91)  d(z)  =  dQ  +  dxz  +  ...  +  d^z”,  iJO 

of  degree  m  >  n.  The  following  theorem  provides  several 
equivalent  characterizations  of  d(z)  and  shows  that 
z^,...,  z^  are  to  be  found  among  its  m  roots. 

THEOREM  5 .  Let  y  =  (v^ . y^ . y ? . . .  )  nave  the  repre¬ 

sentation  (3)  with  a^,...,  a^  being  nonzero  and  with  z^,..., 

z  being  distinct.  Let  coJ  c,,...,  c  be  the  coefficients 
n  U  1  n 

of  the  characteristic  polynomial  c(z)  of  (6)  naving 

T 

z.,...,  z  as  its  roots.  Let  d  =  (d~,d,,...,  d  )  with 
In  Ulm 

m=n+r  and  r  >  0,  and  let  d(z)  be  the  corresponding 
polynomial  (  9 l).  Then  the  following  are  equivalent: 

(i)  d  is  a  (right)  null  vector  for  the  data  matrix 

>0  •  •  •  'm 

(92)  =  *1  y2  •  •  •  "m+l 

•  •  • 

■  *  • 

(  i  i  )  d  (E)  y  =  O, 
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(iv)  d(z)  has  the  factorization 

d(z)  =  (bQ  +  b^z+  ...  +  brzr)  •  c(z) 

for  some  choice  of  the  coefficients  b^jb^,...,  br  , 

T 

(v)  d  lies  in  the  row  space  of  the  (r+l)  x  (m+l)  matrix 


(9  4)  dT  =  bT  •  C 


for  some  choice  of  b1  =  (bQ,  b^...,  br ) 
Proof.  From  the  identity 


tl 


d(E)y  =  dQy  +  d}Ey  +  . . .  +  dmEmy  =  ^ 


we  see  that  (i)  and  (ii)  are  equivalent.  Next,  since  the 

2 

power  sequences  (l,z  ,z  ,...  ),  v  =  1,...,  n  are  linearly 

independent,  and  since 

d(E)yk  =  d(E)  £av*vk  =  £  [avd(zy)  ]zvk,  k:--0, 1 , 2,  .  .  . 


v=l 


V=1 


we  see  that  (ii)  holds  if  and  only  if  avd(zy)=  0  for  each 


5  3 


v  =  1,...,  n  which  is  equivalent  to  (iii)  and  thus  to  (iv) 
since  the  a^'s  are  nonzero.  Finally,  by  equating  the 
coefficients  of  like  powers  of  z  in  the  identity 

d_.+d,  z+.  .  . +d  zm  =  (b-.+b.  z+.  .  .  +b  zr )  (cn+c.  z+.  .  .  +c  zn  ) 
m  Ul  rui  n 

=  bQc ( z )  +  b1zc(z)+  ...  +  brzrc(z) 
we  see  that  (iv)  is  equivalent  to  (v).  Q 

NOTE.  Since  the  conditions  (iii),  (iv),  and  (v)  above 
do  not  depend  upon  any  particular  choice  of  y,  these 
equivalent  conditions  imply  that  (i)  and  ( i i )  must  both 
hold  for  every  possible  nondegenerate  choice  of  y  having 
the  form  (3)  with  the  same  zv's. 

There  are  several  ways  which  can  be  used  to  numerically 
generate  a  family  of  vectors 

(95)  d.T  =  (diQ,  d.m),  i=0,l,...,  r 

of  the  form  characterized  in  Theorem  5.  We  have  already 
observed  that  if  y  is  a  nondegenerate  transient  sequence 
of  the  form  (3)  and  we  perform  a  singular  value  decom¬ 
position  of  the  data  matrix  id  with  m  =  n+r,  then  r  +  1 
of  the  singular  values  will  be  zero  and  the  corresponding 
right  singular  vectors  will  be  null  vectors  of  and 

thus  serve  as  the  ' s  .  In  some  cases,  it  is  convenient 
to  observe  a  number  of  transient  sequences  7t 
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which  result  from  different  initial  excitations  of  the  same 


l 


f 


underlying  physical  system.  Ve  can  then  numerically 
generate  a  null  vector  dL  for  the  data  matrix  of  (  92) 

constructed  using  the  i-th  transient  sequence  y^ , 
i=0,l,...,  r.  The  same  procedure  can  be  used  to  generate 
the  dL  ' s  when  only  one  nondegenerate  transient  y  is  known 
provided  that  we  first  generate  suitable  auxiliary 
sequences  y^  =  P^y,  i=0,l,...,  r  from  y  by  applying 
preprocessing  operators  Pq,  P^,...,  t"1  .  Indeed,  in 
principle,  we  can  map  any  given  nondegenerate  transient 
y  having  component  (3)  into  an  arbitrary  transient  u 
having  components 

n  k 

uk  =  E  av*  zv  »  k=0, 1,2,... 
v=l 

(for  the  same  underlying  system)  by  using  any  preprocessing 
operator  P  =  y(E)  constructed  from  a  function  Y  (z)  which 
is  analytic  on  |z|  <  1  and  which  interpolates  the  points 

Y(zv)  =  av*/av  »  v  =  n  • 

Ve  would  like  to  have  some  numerical  procedure  for 
obtaining  the  system  poles  z^,...,  z^  from  such  a  collection 
of  vectors  (  95).  In  principle,  we  could  simply  factor 
each  of  the  polynomials 

(  96)  d.(z)  =  diQ  +  d^z  +  ...  +  dimzm  ,  i=0,l,...,  r 
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9  •••  i 


and  choose  as  z, 


z  the  roots  which  are  common  to  all 
n 


of  them  (assuming  that  the  dL  ' s  have  been  suitably 
restricted  so  as  to  rule  out  the  possibility  of  an  extraneous 
zero  held  in  common  by  each  of  the  polynomials  (96).)  In 
practice,  however,  the  roots  are  subject  to  slight  pertur¬ 
bations  due  to  noise,  computer  roundoff,  etc.,  and  instead 
of  finding  n  roots  which  are  held  in  common  by  the  poly¬ 
nomials  (  96),  we  obtain  n  clusters  of  roots  (near  z^,..., 
z  )  which  must  be  suitably  averaged.  For  this  reason,  we 
would  like  to  have  some  way  to  process  the  vectors 
dQ,  d^,  ...,  dr  so  as  to  obtain  directly  a  good  estimate 
of  the  coefficients  c  =  (Cq,  c^,...,  c  )  of  the  character¬ 
istic  polynomial  (6)  having  z^,...,  z n  as  roots.  One  very 
nice  scheme  for  doing  this  has  been  published  recently  by 
Henderson  [4],  and  we  shall  now  expand  upon  his  work. 

Suppose  then  that  we  have  been  given  a  collection 
(  95)  with  z^,...,  being  common  roots  of  each  of  the 
corresponding  polynomials  (  96).  By  using  Theorem  5-v  on 
a  row-by-row  basis  we  see  that  the  matrix 


«  doo 

doi 

d0m 

dll 

...  d. 

lm 

j  (  97)  D  = 

1 

.  dr0 

dr  1 

.  .  .  d 

rm 

T  T  T 

(having  dQ  ,  d^  ,...,  d^  as  rows)  can  be  factored  in  tne  form 


(98)  D  =  B  C 


where  B  is  an  (r+l)x  (r+l)square  matrix 

b00  b01  ‘  ‘  *  b0r 

b10  bll  *  ‘  *  blr 

(  99)  B  = 

b  „  b  .  .  .  .  b 
rO  rl  rr 

(having  rows  which  correspond  to  the  polynomial  factors 
of  Theorem  5-iv)  and  where  C  is  the  (r+l)  x  (m+l)  matrix 
(  93).  Our  goal  is  to  extract  the  parameters  Cq,c^,...,  c^ 
of  C  directly  from  the  matrix  D. 

Henderson's  scheme  for  finding  the  c^'s  begins  vith 
the  use  of  Gaussian  elimination  with  partial  pivoting  to 
systematically  zero  out  the  elements  of  the  matrix  D  having 
indices  i,  j  with  i  >  j,  i.e.,  those  elements  whieh  lie 
below  the  principle  or  left  diagonal  drawn  through  d^, 
d.^,  ...  .  This  process  replaces  D  by  the  matrix  LPD  where 

P  is  obtained  by  suitably  permuting  tne  rows  of  tne 
(r+l)  x  (r+l)  identity  matrix  and  where  L  is  an  (r+l)  x 
(r+l)  lower  triangular  matrix  having  l's  along  its  diagonal, 
cf.  [8,  Chapter  l].  The  application  of  P  serves  to  permute 
the  rows  of  D  and  the  subsequent  application  of  L  then  serve 
to  carry  out  the  elementary  row  operations  which  introduce 
the  desired  zero  structure.  Henderson's  scheme  then 
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continues  by  using  additional  row  operations  to  zero  out 
those  elements  of  LPD  having  indices  i,  j  with  j  >  i+m-r, 
i.e.,  those  elements  which  lie  above  the  right  diagonal 
drawn  up  through  d  ,  d  .  .  ,  ...  .  In  this  way  the 

matrix  LPD  is  replaced  by  the  matrix  ULPD  where  U  is  an 
(r+l)  x  (r+l)  upper  triangular  matrix  having  l's  along  its 
diagonal.  The  remaining  matrix  is  then  a  row  scaled 
version  of  C ,  i.e., 

(100)  ULPD  =  SC 
where 

(101)  S 

The  i-th  row  of  ULPD  thus  contains  s.c_,  s.c,,...,  sc 
in  columns  j=  i,  i+1,...,  i+n,  respectively,  so  after  a 
suitable  normalization  (or  averaging  process  cf.  [4, 
p.  986])  c0,  c^,...,  c^  are  obtained.  Henderson  has  shown 
that  this  procedure  will  always  work  when  D  has  full  rank. 

A  slightly  more  general  (necessary  and  sufficient)  condition 
for  the  success  of  this  scheme  is  given  in 

THEOREM  6 .  Let  the  matrix  D  of  (  97)  have  rank 
r+l-p  where  0  <  p  <  r  ,  and  assume  that  »D  has  the  factori¬ 
zation  (  93)  with  B  as  in  (  99)  and  with  C  as  in  (9  3)  with 
Cn^°  *  Then  there  e^ist  (r+l)  x  (r+l)  matrices  P,L,U,S  with 
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P  being  a  permutation  matrix,  with  L  being  a  lower  trian¬ 
gular  matrix  with  unit  diagonal,  with  U  being  an  upper 
triangular  matrix  with  unit  diagonal,  and  with  S  being  a 
diagonal  matrix  such  that  (100)  holds  if  and  only  if 
exactly  p  columns  of  B  vanish  identically. 

Proof :  Suppose  first  that  D,  C,  P,  L,  U,  S  are 

related  as  in  (100).  By  using  (  98)  in  (100)  we  see  that 

(ULPB  -  S)  C  =  0 

and  since  the  last  r+1  columns  of  C  are  linearly  independent 
this  implies  that 

(102)  ULPB  =  S  . 

By  hypothesis,  the  matrices  U,  L,  P  are  nonsingular  so  by 
using  (  98)  and  (102)  we  see  that  D,  B  and  S  must  have 
exactly  the  same  rank  r+l-p  and  that  B,S  have  the  same 
null  space.  In  particular,  exactly  p  of  the  diagonal 
elements  s^  of  S  vanish,  the  corresponding  p  columns  of 
B  must  vanish,  and  (since  the  rank  of  B  is  r+l-p)  no  other 
columns  of  B  can  vanish. 

Conversely,  assume  that  exactly  p  columns  of  B  vanish 
(with  r+l-p  being  the  rank  of  both  D  and  B  ).  By  using 
elementary  row  operations  we  can  reduce  B  to  an  upper 
triangular  matrix  having  exactly  p  zeros  along  the  principle 
diagonal,  i.e.,  we  can  find  P,L  such  that  LPB  is  an  upper 
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triangular  matrix  with  exactly  p  zero  elements  along  the 
diagonal.  Since  each  zero  column  of  B  is  also  a  zero 
column  of  (LP)  •  B,  the  p  columns  of  LPB  which  contain  the 
p  zero  pivots  must  vanish.  This  being  the  case,  elementary 
row  operations  can  be  used  to  reduce  LPB  to  a  diagonal 
matrix,  i.e.,  we  can  find  U,S  such  that  (102)  holds.  Upon 
multiplying  (102)  on  the  right  by  C  and  using  (  98)  we 
then  obtain  (100).  {] 

NOTE;  When  D  satisfies  the  conditions  of  the  theorem 
the  matrices  LPD  and  ULPD  =  SC  will  be  obtained  naturally 
during  the  two  stage  elimination  process.  Indeed,  suppose 
that  D  can  be  factored  in  the  form  (  98)  and  that  Cq^  0. 

(If  Cq  =  0  we  see  from  (103)  and  (  98)  that  the  first 

columns  of  C  and  D  both  vanish  and  that  ve  could  replace 
our  problem  by  one  corresponding  to  a  smaller  value  of  r. ) 
Since  D  has  rank  r+l-p  ,  when  we  use  elementary  row 
operations  to  reduce  D  to  upper  echelon  form  we  will  end 

up  with  exactly  p  zero  rows.  From  (9  8)  we  see  that  each 

row  of  D  is  some  linear  combination  of  the  row’s  of  C,  and 
since  c^O  we  see  that  each  of  the  r+l-p  nonzero  rows  of 
our  upper  echelon  matrix  must  have  at  least  one  nonzero 
element  among  its  first  r+1  components.  This  being  the 
case  we  can  rearrange  the  p  zero  rows  so  as  to  obtain  the 
matrix  LPD  which  has  zeros  below  the  principle  diagonal  and 
exactly  p  zeros  (which  lie  in  the  p  zero  row’s)  along  the 
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principle  diagonal.  Thus  we  see  that  the  first  stage  of  the 
elimination  process  can  be  carried  out  on  any  matrix  D 
having  the  factorization  (  98). 

It  is  the  second  stage  of  the  elimination  process  that 
is  dependent  upon  the  auxiliary  hypothesis  that  B  has  exactly 
p  zero  columns.  Indeed,  if  the  (lower  right)  r,m  -  element 
of  LPD  is  nonzero,  then  by  subtracting  suitable  multiples 
of  thi s  row  from  the  previous  ones  we  can  zero  out  the 
upper  most  r  elements  of  the  last  column.  If  this  r,m  - 
element  is  zero,  however,  we  can  proceed  if  and  only  if  the 
whole  last  column  of  LPD  is  already  filled  with  zeros.  Since 

LPD  =  (LPB)  •  C 

we  see  that  the  last  column  of  LPD  will  vanish  if  and  only 
if  the  last  column  of  C  is  in  the  null  space  of  LPB  and 
thus  in  the  null  space  of  B,  i.e.,  if  and  only  if  the  last 
column  of  B  vanishes.  Analogous  considerations  apply  at 
subsequent  stages  of  the  back  elimination  process. 

NOTE.  If  the  matrix  D  has  full  rank  (as  is  often  the 
case  in  practice)  then  the  square  matrix  B  from  (  98)  must 
also  have  full  rank  r+1  so  that  p=0  and  none  of  the  columns 
of  B  vanish.  Theorem  6  then  guarantees  that  D  has  the 
factorization  (lOO).  Moreover,  in  this  case  S  also  has 
full  rank  so  that  U,L,S  are  uniquely  determined  by  D  and 
the  pivoting  strategy  P. 
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The  following  example  has  been  constructed  to  show 
that  there  are  cases  where  the  two  stage  elimination  process 
of  Theorem  6  fails  even  when  the  roots  z^,...,  z are 
uniquely  determined  by  d^,  d^,...,  df  .  Indeed,  let 

d^(z)  =  l+2z+z^+2z^  =  2(z-  l/2) (z-i ) (z+i ) 

d^(z)  =  z^+  z^  =  ( z-0) ( z-O) ( z-i ) ( z+i ) 

d2(z)  =  l+2z+2z^+2z^+z4  =  ( z-1 ) ( z-1 ) ( z-i ) ( z+i ) 

so  that  d^Cz),  d^(z),  ^2(z)  have  only  the  roots  z  =  +  i  in 
common.  In  this  case 


1  2  0 

10  10  0 

B  = 

0  0  1 

C  = 

0  10  10 

.12  1 

* 

0  0  10  1 

•* 

and 


D  =  BC 


12  12  0' 
0  0  10  1 
1  2  2  2  1 


The  matrices  B,D  both  have  rank  2  so  that  p  =1,  but  no 
column  of  B  vanishes.  Upon  carrying  out  the  forward 
elimination  process  on  D  we  find 


1 

2 

1 

2 

o' 

'  1 

2 

1 

0 

i— 

O1 

1 

2 

1 

2 
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0 

0 

1 

0 

1 

-> 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 

1 

2 

2 

2 

1 

0 

0 

1 

0 

1 , 

0 

0 

1 

0 

1 

The  back  elimination  fails,  however,  at  the  point  where 
we  try  to  annihilate  the  2  in  the  s ec ond-to- 1  a s t  column 
of  the  first  row. 
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8.  CONCLUSIONS 


Our  approach  to  the  problem  of  transient  analysis, 
i.e.,  the  problem  of  extracting  the  system  poles  z^...,  z^ 
from  a  sequence  of  samples  y  =  (yQ>  y^,  y2»...)  involves  a 
three  stage  process.  Ve  first  assemble  a  data  matrix  by 
applying  certain  sequence-to-sequenee  mappings  to  y 
within  the  conceptual  framework  of  Theorem  3  or  Theorem  4. 

Ve  then  compute  a  null  vector  for  this  data  matrix  by  using 
an  eigenvalue  analysis  or  a  singular  value  decomposition 
(with  the  former  being  less  costly  of  computer  storage  and 
the  latter  being  somewhat  better  conditioned.)  This  null 
vector  then  yields  a  characteristic  function  having  z^ , .  .  .  , 
zn  as  roots.  Our  scheme  is  a  conceptually  simple  one  which 
admits  significant  new  generalizations  (such  as  those  of 
Examples  4  and  7  in  Section  6),  and  it  places  the  existing 
algorithms  within  a  common  mathematical  framework. 

The  sequenc e-to-sequenc e  mappings  which  lie  at  the 
heart  of  our  analysis  can  be  effected  quite  simply  on  a 
digital  computer.  Ve  have  analyzed  the  noise  suppressing 
properties  of  such  ma.ppings  and  identified  a  fundamental 
trade  off  between  the  size  and  shape  of  the  region  of 
amplification  and  the  signal-to-noi se  ratio  improvement 
which  can  be  achieved.  The  analysis  clearly  shows  why  Jain's 
method  is  successful  in  filtering  noise  from  the  low  order 
poles  when  high  sampling  rates  are  used  (i.e.,  when  the  first 
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few  zv ' s  are  close  to  l)  and  points  the  way  to  other  pre¬ 
processing  schemes  which  will  facilitate  the  computation 
of  the  higher  order  poles.  Further  noise  reduction  could 
be  achieved  by  using  an  adaptive  scheme  which  first  estimates 
the  system  poles  and  then  carefully  computes  them  using 
preprocessing  operators  which  are  optimal  with  respect  to 
the  expected  pole  pattern. 

Finally,  the  problem  of  estimating  the  system  order  n 

and  assessing  the  accuracy  of  the  computed  poles  z^,...,  z^ 

can  also  be  solved  by  using  pole  preserving  mappings.  In 

the  absence  of  noise  the  given  sequence  y  and  the  auxiliary 

sequences  yQ  =  Y0(E)y,  =  Y1(E)y, ...,  y r  =  Yr(E)y  will  all 

have  exactly  the  same  poles  z^,...,  with  n  unknown.  Ve 

can  perform  separate  computations  of  n  and  the  zv's  using 

each  of  the  y^'s,  or  we  can  compute  approximate  null 

vectors  d^,  d^,...,  dv  for  data  matrices  constructed  using 

Tr\t  3Ti  *  •  .  .  »  y.  and  then  extract  n  and  the  z  '  s  from  the 
u  “i  r  v 

d^'s  by  using  Henderson's  method. 

Many  of  these  concepts  have  been  tested  by  performing 
the  related  computations  on  simple  examples,  and  the  results 
have  been  most  encouraging.  A  good  deal  more  work  will  be 
required,  however,  in  order  to  perfect  these  ideas  and  to 
incorporate  them  into  efficient  production  codes. 
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9.  APPENDIX  -  THE  DAMPED  VIBRATING  STRING 

Problems  in  transient  analysis  arise  when  a  physical 
system  reverberates  in  response  to  some  initial  excitation. 

In  principle,  such  a  phenomenon  can  usually  be  modeled  by 
solving  a  certain  boundary  value  problem  (which  characterizes 
the  system)  subject  to  initial  conditions  which  depend  upon 
the  form  of  the  initial  stimulus.  In  practice,  such  models 
can  be  analyzed  in  detail  only  in  extremely  simple  situations 
where  there  is  unusual  symmetry  or  low  dimensionality.  The 
damped  vibrating  string  provides  us  with  a  simple  pnysical 
phenomenon  which  nicely  serves  to  illustrate  the  way  a 
problem  in  transient  analysis  arises  and  the  inherent 
difficulties  associated  with  its  solution. 

Let  u(x,t)  give  the  (one  dimensional)  displacement 
from  the  equilibrium  position  of  the  string  at  coordinate 
x,  0  <  x  <  L,  at  time  t  >  0.  The  motion  of  the  string  is 
governed  by  the  partial  differential  equation 

(103)  t  uxx(x,t)  =  p  utt(x,t)  +  a  ut(x,t),  0<x<L,  t>0 

where  t  is  the  tension  of  the  string,  p  is  the  linear  mass 
density,  and  *  is  the  damping  coefficient.  Ve  shall  assume 
the  endpoints  of  the  string  are  fixed  by  the  boundary 
conditions 

(104)  u(0,  t )  =  u ( L,  t )  =0,  t  >  0, 
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and  that  some  external  stimulus  has  subjected  the  string  to 
the  initial  displacement  and  velocity 

(105)  u(x'°)  =  u0<x>'  0  <  x  <  L 

ut (x» 0)  =  vQ(x),  0  <  x  <  L 

at  time  t=0.  Upon  separating  variables  ve  find  that  any 
solution  of  (103)  -  (104)  has  the  form 

(106)  /  .  i  x^3  <  k  (-a+it u  )t,-7  (-a-iwm)t)  .  ,  „  \ 

u(x,t)=  2^  |  A^e'  nr  +  Ame  m  jsin(mrtx/L) 

m=l 

2 

where  i  =  -1  and  where 
a  =  */2p 

l.2p  4p2  ) 

are  given  in  terms  of  the  physical  parameters  t ,  p,  k,  L 
which  characterize  the  system.  Vhen  (106)  is  subjected  to 
the  initial  conditions  (105)  we  find  that  the  complex 
coefficients  A-^,  A.2>  ...  are  given  by  the  integrals 
L 

Am"^U,mL^~1  /  1  %)U0 ( x )  ~  i  t  auQ  ( x )  +vo  ^ x  ^  sin(mT1x/L)<lx,m=l,2,  •  ’  • 
x=0 

and  thus  depend  on  the  choice  of  Uq(x),  Vq(x). 

Suppose  now  that 

Y(t)  =  u(xQ, t)  j  t  >  0 
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as  the  resulting  da  spl  a  e  ament  of  1 1  j  e  string  at  some  fixed 
point  Xq,  0  <  Xq  <  L.  Using  (106)  ve  see  that  the  transiem 
Y(t)  can  be  written  in  the  form 

T  ( t  )=  [Amsin(mnx0/L)  ]e^-a+ia,m^t-‘-[Amsin(mnx0/L)  ]e^-a_iu,m^1  ] 

m=l 

or  equivalently  in  the  form 
cc- 

T(t)  =  £  av  eSvt 

v=l 

where 

ay  =  Amsin  (mnxQ/L)  and  sy=  -a+iu)m  if  v=2m-l ,  m=l,2,... 

a>,  =  A  sin  (m-x./L)  and  s  =  -c-ia  if  v=2m.  m=1.2,...  . 

J  iTi  u  *  m 

The  poles  sv  are  more  or  less  regularly  spaced  along  the 
vertical  line  Rez  =  -a  in  the  left  half  plane.,  cf.  Fig.  1,  and 
they  depend  only  on  the  physical  parameters  ' . p . * , L  of 
the  system.  On  the  other  hand,  the  coefficients  av  depend 
on  the  point  of  observation  Xq  and  the  initial  excitation 
(i.e.,  Uq(x)  and  v^(x)  .)  The  problem  of  transient  analysis 
is  then  to  recover  the  sv's  from  certain  samples 

yk  =  I(kT),  k=0, 1,2,... 

of  the  signal  Y(t). 

The  intrinsic  difficulty  of  the  problem  is  now  apparent. 


At  best  we  can  hope  to  find  the  first  few  of  the  infinitely 


Moreover,  if  the  initial  stimulation 


many  system  poles,  sv  . 
of  the  string  fails  to  excite  a  given  mode,  if  the  point  of 
observation  Xq  happens  to  lie  at  a  node  of  that  mode,  or 
if  we  unwittingly  choose  too  small  a  sampling  rate,  then 
the  mode  will  be  weakly  represented  (if  at  all)  in  the 
sequence  of  samples  y  =  (y^,  y^>  y2>...  )  an<i  we  viH  fail 
to  find  the  corresponding  sv  .  Nevertheless,  in  practice 
we  find  that  if  we  use  a  reasonable  sampling  rate  then  it 
is  possible  to  extract  at  least  the  first  few  sv's  for 
"almost"  all  choices  of  x^,  Uq(x),  and  v^(x)  •  Analogous 
considerations  apply  when  we  use  transient  analysis  to 
study  more  complex  physical  systems  which  cannot  be 
subjected  to  such  a  detailed  analysis. 
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(a)Small  T 


(b)  Intermediate  T 


(c)  Large  T 


Figure  1.  The  transformation  z=exp(sT)  of  (4)  maps  the  sv's 
from  the  left  half  s-plane  onto  the  zv ' s  within  the  unit 
circle  of  the  z-plane,  with  T  >  0  being  the  sampling  interval 

(a)  When  T  is  too  small,  the  zv's  are  tightly  clustered 
near  z=l. 

(b)  When  T  is  correctly  chosen,  the  first  few  zv's  are 
nicely  separated. 

(c)  When  T  is  too  large,  tne  z^'s  are  hurried  deep  within 
the  unit  circle. 


Figure  2.  Contour  plots  of  the  SNR  amplification  factor 

2  2  2 
g(z)  =  |  1+z+z  |  /3  for  the  case  where  y(z)  =  1  +  z+z 

Tne  region  of  amplification  CL  is  the  set  of  points  z 

within  the  unit  circle  for  wnich  g(z)  >  1. 


Figure  3.  The  region  of  amplification  CL  for  the  case  where 
Y(z)=l+z+  ...  +z^  ^  and  N=2,3,4,10.  As  N  increases,  the 
maximum  SNR  amplification  g(l)=N  also  increases  but  CL 
shrinks  in  size. 


Figure  4.  The  region  of  amplification  for  the  case  where 
Y(z)  =  l+z/R  +  z  /IT  +  ...  and  R=2,  1.4,  1.1,  1.01.  As 
R  ^  1,  the  maximum  SNR  amplification  g( 1 )=(R+1 )/(R-l ) 
increases  but  CL  snrinas  in  size. 


Figure  5.  The  region  of  amplification  for  the  case  where 
y(z)  is  given  by  (74)  with  p=2  and  with  the  four  poles 
(R^j  =  as  shown.  In  this  case  CL  has  the 

approximate  form  (71)  and  the  corresponding  preprocessing 
scheme  is  well  suited  for  pole  patterns  of  the  form  shown 
in  Fig.  1  (b). 
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